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ABSTRACT 

Inclusion of moist physics in the linearized version of a weather forecast model is beneficial in terms of 
variational data assimilation. Further, it improves the capability of important tools, such as adjoint-based 
observation impacts and sensitivity studies. A linearized version of the relaxed Arakawa-Schubert (RAS) 
convection scheme has been developed and tested in NASA’s Goddard Earth Observing System data as- 
similation tools. A previous study of the RAS scheme showed it to exhibit reasonable linearity and stability. 
This motivates the development of a linearization of a near-exact version of the RAS scheme. Linearized 
large-scale condensation is included through simple conversion of supersaturation into precipitation. The 
linearization of moist physics is validated against the full nonlinear model for 6- and 24-h intervals, relevant to 
variational data assimilation and observation impacts, respectively. For a small number of profiles, sudden 
large growth in the perturbation trajectory is encountered. Efficient filtering of these profiles is achieved by 
diagnosis of steep gradients in a reduced version of the operator of the tangent linear model. With filtering 
turned on, the inclusion of linearized moist physics increases the correlation between the nonlinear pertur- 
bation trajectory and the linear approximation of the perturbation trajectory. A month-long observation 
impact experiment is performed and the effect of including moist physics on the impacts is discussed. Impacts 
from moist-sensitive instruments and channels are increased. The effect of including moist physics is exam- 
ined for adjoint sensitivity studies. A case study examining an intensifying Northern Hemisphere Atlantic 
storm is presented. The results show a significant sensitivity with respect to moisture. 


1. Introduction 

The techniques and methods used in incremental var- 
iational data assimilation are largely based on the as- 
sumption that the underlying behavior of the system is 
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linear or close to linear (Courtier et al. 1994). To apply 
the techniques, such as minimization of a cost function in 
four-dimensional variational data assimilation (4DVAR), 
a linearized version of the model is used. 

The linear model can be thought of as providing an 
approximation to the trajectory of nonlinear perturba- 
tions. For example, it may be used to estimate how much 
an initial error (perturbation from the truth) would grow 
over time. If the full model is in fact linear, then the 
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approximation would give the exact nonlinear pertur- 
bation trajectory. 

In practice, how well a linearized model will perform 
for data assimilation applications depends on three main 
issues. First, how much detail is lacking from the linear 
model, for example physics that is not included; second, 
how close to linear the components of the nonlinear 
model are; and, third, how many approximations are 
made in the numerics — for example, use of multiple 
outer loops without higher than first-order lineariza- 
tion techniques (Tremolet 2008). 

For the large-scale dynamics of the atmosphere a lin- 
earization produces a very good approximation of the 
quasi-linear perturbations. Indeed, for a primitive equa- 
tion model, with only very simple physical parameteri- 
zations, a correlation of 0.9 or more between the linear 
and nonlinear perturbation trajectories is possible, even 
for 72-h integrations (Errico et al. 1993). Atmospheric 
motions that are subgrid scale, sometimes referred to 
as “the physics,” include processes such as turbulence, 
convection, precipitation, gravity wave drag, and ra- 
diation. The physical parameterizations that are used 
to model these processes can be very nonlinear and 
contain discontinuities. Including these processes in the 
linear model in a way that produces accurate represen- 
tation of the actual perturbation trajectory requires much 
care. 

The National Aeronautics and Space Administra- 
tion’s (NASA) Global Modeling and Assimilation Of- 
fice (GMAO) is currently developing a linearization of 
the Goddard Earth Observing System (GEOS-5) atmo- 
spheric global circulation model (AGCM). The linear 
model is based on the new cubed-sphere dynamical core 
(Putman 2007) and replaces the old latitude-longitude- 
based linearized model. So far, the adjoint and tangent 
linear versions of the dynamical core have been de- 
veloped, along with a simple vertical diffusion scheme 
and boundary layer. The new linearized model pro- 
vides an essential upgrade to the adjoint-based obser- 
vation impact tool, used for daily monitoring. It has 
also been designed to provide a 4DVAR-capable as- 
similation system. 

The analysis produced by a 4DVAR system can be 
significantly improved if the assimilation of observations 
affected by moist processes such as clouds and pre- 
cipitation is possible (Amerault et al. 2008; Errico et al. 
2007; Errico and Raeder 1999; Janiskova et al. 1999; 
Lopez and Moreau 2005; Mahfouf and Rabier 2000; 
Stiller and Ballard 2009; Stiller 2009; Tompkins and 
Janiskova 2004). This requires the inclusion of moist 
physics in the adjoint and tangent linear models. In- 
cluding moist physics will also improve the estimation 
of the impacts coming from instruments such as the 


Atmospheric Infrared Sounder (AIRS), High Resolution 
Infrared Radiation Sounder (HIRS), and Microwave 
Humidity Sounder (MHS), as well as feature-tracking 
satellites, all of which have a sensitivity to moist pro- 
cesses. Further applications that can benefit from having 
moist physics in the linearization are adjoint-based sen- 
sitivities (Jung and Kim 2009) and singular vector cal- 
culations (Ehrendorfer and Errico 1995). Additionally, 
three-dimensional variational data assimilation (3DVAR), 
4DVAR, and many operational ensemble methods 
make use of the linear model in the observation operator. 
Including accurate moist physics in the linear model is 
essential when assimilating moisture-affected observa- 
tions, such as those that will be available from the up- 
coming Global Precipitation Satellite (Hou et al. 2008). 

Convection in the nonlinear GEOS-5 AGCM is mod- 
eled using the relaxed Arakawa-Schubert (RAS) con- 
vection scheme (Arakawa and Schubert 1974; Moorthi 
and Suarez 1992). Large-scale condensation is modeled 
using the scheme developed by Bacmeister et al. (2006). 
The convection is computed prior to the large-scale 
condensation. Holdaway and Errico (2013) examined 
linearity and stability in the RAS convection scheme by 
studying Jacobian sensitivities. They found the scheme 
to generally exhibit good linearity and stability prop- 
erties. Based on the findings of that work, a few simple 
modifications are applied to the RAS scheme and then 
an exact linearization is developed. The Bacmeister 
et al. (2006) scheme is rather more complex than the 
RAS scheme. It contains strong nonlinearities and re- 
lies on a large number of inputs that are not readily 
available in the linearized model. Rather than attempt- 
ing to implement an exact linearization of this scheme, 
a reduced large-scale condensation scheme is imple- 
mented, one that simply converts supersaturation to 
precipitation and warming (Errico et al. 1994). 

Constructing the linearized moist physics and exam- 
ining its behavior against the full nonlinear system is 
very useful for understanding how the sensitivities in 
the moist schemes behave. This can assist not only in 
improving the analysis and associated data assimilation 
tools but also in developing a general understanding 
of how the system responds to sensitivities in the moist- 
physics schemes. This is important for anyone devel- 
oping data assimilation systems. Even ensemble Kalman 
filter type methods, which may not directly require the 
linear model, do rely on certain assumptions about line- 
arity in the system. The developers of the full nonlinear 
moist-physics schemes can also benefit from the iden- 
tification of sensitivities that are not realistic and thus 
where the scheme may require improvement. 

The development of the linearized moist-physics 
schemes is outlined in section 2. Validation of the linear 
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approximation for the moist physics is examined in 
section 3. The effect that moist physics has on the ob- 
servation impacts is discussed in section 4. In section 5 
a case study of an intensifying Atlantic storm is used to 
discuss the impact of linearized moist physics on ad- 
joint sensitivity studies. Section 6 discusses the role of 
moisture in the norm. Section 7 offers some concluding 
remarks. 

2. Development 

In this section the general approach to modeling the 
linearized moist physics is outlined. Modern numerical 
schemes, including the RAS scheme, are complex it- 
erative procedures so the full linearization is not pre- 
sented. Instead, some key components of the schemes 
being used are outlined, any simplifications added to 
the schemes are listed, and any linearization issues are 
addressed. 

The prognostic variables that will be used in the moist- 
physics schemes in the GEOS-5 linear model are the 
zonal and meridional wind speeds u and v (m s -1 ), the 
potential temperature 6 (K), the specific humidity 
q (kg kg” 1 ), and the surface pressure p s (Pa). The pro- 
files of temperature and specific humidity, along with the 
surface pressure, are used to parameterize the moist 
processes that are occurring. For the purposes of this 
work the surface pressure is not altered by the moist 
processes. The horizontal wind speeds can be altered, 
for example through vertical transport of momentum by 
convection, but are not explicitly used in determining 
the moist processes that are occurring. In the following 
text when linearizing variables, a superscript prime (') is 
used to denote perturbation parts and a superscript r is 
used for the reference parts (e.g., 6 = 0 <r> + 6'). 

a. Convection 

The relaxed Arakawa-Schubert scheme was first in- 
troduced by Moorthi and Suarez (1992). The scheme of- 
fers a number of simplifications on the original Arakawa 
and Schubert (1974) scheme to make it suitable for use 
in general circulation models. The main simplification 
is that the system is relaxed toward equilibrium at each 
step, as opposed to being iterated to a fully equilibrated 
state every time step. 

Convection in the atmosphere is an inherently non- 
linear process. Mechanisms that need to be represented 
include latent heating, deposition, fast updrafts and 
downdrafts, and reevaporation. In addition to this the 
schemes used to represent convection often make use 
of discontinuous modeling and employ artificial pro- 
cesses that increase nonlinearity. A simplified expla- 
nation of the RAS scheme is given here to demonstrate 


the linearization limitations. See Moorthi and Suarez 
(1992) for a full derivation and explanation of the scheme. 

The RAS scheme updates a single column of the at- 
mosphere at one time and considers an ensemble of 
cloud depths within the column. Each cloud depth is 
assumed to have its base at the same model level, cur- 
rently given by the top of the boundary layer. Different 
clouds have different detrainment levels, starting from 
the base of the cloud up to a maximum (around 30 hPa 
in the current formulation). The model has 72 levels in 
total and the lid is at 0.01 hPa. Clouds in the ensemble 
are characterized by an entrainment parameter A. The 
normalized mass flux for each cloud depth is linear with 
height. Central to the algorithm is the calculation of a 
cloud work function for each cloud depth. The RAS 
scheme models all depths of convection using the same 
algorithm. Shallow convection just has fewer members 
since clouds detraining at higher levels are not found. 

Before the effects of individual cloud depths are ap- 
plied to the atmospheric profile, the RAS scheme de- 
termines whether convection is occurring (i.e., whether 
cloud is detraining at that model level). This calculation 
is done through a number of conditional statements 
based on the atmospheric profile in that column. In 
all, there are six conditional statements used to check 
whether a cloud depth should be included. The first 
considers whether the relative humidity at the cloud 
base is above a certain threshold. Second, the moist 
static energy at the base layer must be greater than the 
saturation moist static energy at the detraining level. 
Third, the entrainment parameter A [Eq. (A18) in 
Moorthi and Suarez (1992)] must be positive. The en- 
trainment parameter must also not be above a critical 
value, currently 10” 4 kgm” 2 s -1 . The RAS scheme com- 
putes a cloud work function [Eq. (A22) in Moorthi and 
Suarez (1992)], and it has to be above a critical value. 
The critical value itself is modeled discretely based on 
the profile. The final condition is that the rate of change 
of the cloud work function is negative and that the liquid 
water mixing ratio of the detraining air is positive. 

In addition to these explicit conditional statements, 
used to determine the presence of convection, condi- 
tional statements exist in the formulas themselves. 
These discrete steps in the function are evident in the 
algorithm presented in the appendix in Moorthi and 
Suarez (1992) and are represented in the numerics using 
minimum and maximum statements. Discrete modeling 
will result in nonlinearity in the system and inaccuracy 
in a linear version of the scheme. 

When attempting to represent moist physics in the 
linear model, a number of options exist, as outlined 
by Holdaway and Errico (2013). Ideally, an exact lin- 
earization of the schemes would be implemented. 
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However, an exact linearization may quickly diverge 
from the true perturbation trajectory since disconti- 
nuities and nonlinearities in the model are not properly 
taken into account. Alternatively, one could develop 
a new simple scheme to linearize (Lopez 2007). Or one 
could first apply simplifications or smoothing to the 
nonlinear scheme to refine the behavior of the sub- 
sequently derived linear model. All options come with 
challenges and difficulties. 

To estimate the nonlinearity in the RAS scheme, 
Holdaway and Errico (2013) performed an extensive 
study examining Jacobian sensitivities. Despite the po- 
tential nonlinearity in the RAS formulation, they found 
that a large degree of linearity exists, especially for mid- 
to deep convection. In addition the linearized RAS 
scheme was found to be largely stable such that any 
solutions exhibiting growth are relatively well behaved. 

The findings of Holdaway and Errico (2013) motivate 
an approach that involves the exact linearization of 
the RAS scheme, with only minor simplifications to the 
algorithm. These simplifications are related to the way 
certain constants are chosen and have a minimal effect 
on the overall behavior of the scheme. If the highest 
detraining level is less than eight model levels (approx- 
imately 80-100-hPa depth) above the cloud base layer 
and no convective precipitation is occurring, the linear 
RAS scheme is not invoked. The scheme supports 
the carrying of tracers, which can be neglected in the 
linearization. 

A complication when developing a linearization of the 
RAS scheme is due to the use of an ensemble of cloud 
depths, detraining at each model level. As the scheme 
iterates through the cloud depths, the atmospheric pro- 
file is updated. The order is from shallowest cloud to 
deepest. The updated profile is used in the calculation 
for the next cloud depth, and so on. When the adjoint is 
implemented, it will work in reverse order; that is, it will 
start with the deepest cloud in the ensemble and finish 
with the shallowest. To obtain the trajectory used in 
the calculation of the deepest cloud requires iterating 
though each proceeding cloud in the nonlinear model. 
This embedded loop results in either extensive recom- 
putation or large use of memory for the saving of all 
variables for each cloud depth. Both approaches are 
considered here. The optimum approach is found to be 
a combination of recomputation and the saving of vari- 
ables that are computed for every cloud depth. Recording 
the results of all discrete switches avoids unnecessary 
recalculations when iterating the adjoint. 

b. Saturation specific humidity 

Moist-physics schemes rely on the computation of 
the saturation mixing ratio. This is an important quantity 


that describes the mass of water vapor contained in a unit 
mass of saturated air. Mathematically, it is expressed as 


% = e i 


v p(T) 

p~ e 2 v P ( T y 


a) 


where v p (T ) is the saturation vapor pressure and e x = 
0.622 and e 2 = 0.378 are dimensionless constants. 

Saturation vapor pressure is computed using the 
Clausius-Clapeyron equation, a function of tempera- 
ture that involves an exponential. The presence of the 
exponential makes the calculation expensive so in prac- 
tice the formula is replaced with a table lookup. Values 
for v p (T) are interpolated from the table elements: 

” P (T) = v p (TJ+^^[v p (T m+1 )-v p (TJ]. ( 2 ) 


The table elements are denoted with subscript m; T m 
and T m+1 are the nearest values below and above T 
used in obtaining v p from the table. The resolution of 
the table is AT. This piece-wise linear table-lookup 
approach makes the formulation piece-wise linear in T. 
The linearized version of Eq. (2) is 

v P W' = ^ v P V%+ 1 )-^ 7 #)]- < 3 ) 


In the current model the table lookup uses different 
formulations for ice and liquid phases. The table reso- 
lution is AT = 0.1 K and values are computed between 
T = 150.0 and 333.0 K. See Murphy and Koop (2005) 
for a recent review of the methodologies used when 
calculating saturation vapor pressure tables. 

c. Large-scale precipitation 

For the large-scale precipitation a simple scheme that 
precipitates supersaturation, and as used by Errico et al. 
(1994), is applied. The scheme removes supersatura- 
tion at any model level by precipitating out excess 
water while heating the air until the relative humidity 
becomes 1. The amount of latent heating is applied to 
the potential temperature 0. The formulation for large- 
scale nonconvective adjustment is 


7„+i = - a 7 and 


( 4 ) 


0 ,, =0 + 
n + 1 n 


L 


77 C 


A q. 


pm 


( 5 ) 


where subscript n denotes the time step. The moisture- 
dependent specific heat capacity is given by c pm = 
c p { 1 + 0.887(7), and the constant specific heat capacity is 
c p = 1004.49 J kg - 1 K -1 . The latent heat of condensation 
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is L = 2.5104 X 10 6 Jkg -1 . The Exner pressure tt relates 
potential temperature to temperature, T = tt6. 

The profile adjustment is given by 


See Errico et al. (1994) for the derivation of these 
formulas. 

The linearization of Eq. (4) is 


Aq 


(7 - 7 5 ) 
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-1 
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otherwise . 


( 8 ) 


Terms including c:' pm are considered small and are ne- 
glected. The ( dqJdT )' term is included here but is 
sometimes also neglected, as in Errico et al. (1994). 

Note that the formulation above includes a switch, 
when q = q s . This produces a nonlinearity. Away from 
the switch, the functions are linear and will perform well in 
the linearization. Close to the switch, the linear approxi- 
mation may be less accurate. Smoothing could be applied 
near the switch, though this is not considered here. 

The nonconvective precipitation rate is given by 

R b=^Z*qk?k’ ( 9 ) 

where g is the acceleration due to gravity, At is the model 
time step, cr k is the air mass in level k, and A q k is the 
excess moisture in level k. Units are millimeters per 
unit area per second. 

A simplification in the model is that the reevaporation 
of convective and nonconvective precipitation is not 
considered. 

d. Dynamically adjusted trajectory 

It is clear from Eq. (4) that supersaturation caused by 
the dynamics or convection is removed by the scheme. If 
the scheme was run twice in succession, the second run 
would produce no change to the variables. 

An approximation in the system architecture in the 
linearized GEOS-5 AGCM is that the nonlinear tra- 
jectory is not passed between individual components 
of the linear model. For example, as the linearized dy- 
namics is run, the trajectory could be updated and passed 
to the linearized physics components; however, only 
the perturbation trajectory is passed. The nonlinear 
trajectory is read in once per time step (20 min) and all 
components see the same trajectory profile, irrespective 
of the order in which components are called. This 


approximation is suitable when only large-scale dynamics 
and simple boundary layer physics are considered but 
could be problematic once equilibrium-seeking moist- 
physics schemes are included. 

There is a discrepancy between the way that the 
physics and dynamics coupling are handled in the non- 
linear model and the way they are handled in the linear 
model. In the nonlinear model each individual physics 
routine (turbulence, moist physics, gravity wave drag, 
radiation, etc.) produces a tendency for the temperature 
field. These tendencies are then weighted and combined 
at the end of a time step. The specific humidity is coupled 
consecutively; that is, each physics component sees spe- 
cific humidity that has been adjusted by the previous 
component. In the linear model perturbation temperature 
and specific humidity are both handled consecutively. In 
a time step of the tangent linear model the dynamics is 
called and then the moist physics followed by the turbu- 
lence; for the adjoint model, the order is reversed. 

Since the system does not pass trajectories and there is 
a discrepancy in the coupling, a new moist-trajectory 
component is added. The moist components are poten- 
tial temperature and specific humidity, output by the 
nonlinear model just prior to the convection being 
called. This ensures that when the linearized moist 
physics is invoked, it sees unadjusted profiles and pro- 
duces the correct effect. It is sufficient to use the original 
wind and pressure trajectory. The index of the cloud- 
base layer and the highest level of convection are also 
output into the trajectory by the nonlinear model. 
These quantities are required by the RAS scheme and 
are used to determine if convection is deep enough to 
be considered in the linear model. 

e. Energy norm 

Many applications of the linear model require a choice 
of metric. For example, when observation impacts are 
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Fig. 1. Sensitivity initialized at 0000 UTC 18 Mar 2012 and integrated for 6 h using the adjoint 
model. Shown is the sensitivity of dry total energy with respect to the specific humidity dJIdq 
(kgkg -1 ) at 500 hPa. The linearized moist physics is switched on. 


computed, the adjoint is used to determine the reduction 
of the forecast error at the observation time. That error 
must be measured in a certain way. 

A common choice for the metric is the total pertur- 
bation energy, a sum of kinetic energy and a form of 
approximate available potential energy. So, observation 
impacts describe the ability of a certain set of instru- 
ments to reduce the error measured in terms of the total 
energy. The choice of metric is somewhat arbitrary and 
must be carefully chosen with the application in mind. 
For example, Hoover and Morgan (2011) use an ad- 
joint model to examine cyclone track sensitivities. They 
use a vorticity error measure since that is deemed the 
important quantity; using total energy may not describe 
the sensitivities that are most relevant. In another study, 
Errico and Vukicevic (1992) initialize adjoint sensitivity 
integrations with the error in the forecast of surface 
pressure. 

Currently, a dry total energy norm is used in GMAO’s 
operational observation impact tools (Gelaro et al. 2010). 
This is a suitable choice given the dry physics in the cur- 
rent linear model. The interpretation of the total energy 
norm is presented in Errico (2000). In this study the 
use of a moist component in the energy norm metric is 
also considered. 

The formulation that is used for the total energy norm 
follows that given by Ehrendorfer et al. (1999): 


a! 2 + v’ 2 + 


f T 1 + RT n 


1 2 


L 2 

+ E C t (I 
L p 1 o 


/ 2 


( 10 ) 


In Eq. (10), n! and v’ are the perturbed horizontal wind 
components and p' s is the perturbed surface pressure. 


The constants T 0 and p 0 are 270.0 K and 1000.0 hPa, 
respectively, and R = 287.00 J kg -1 K _1 is the gas con- 
stant of dry air. The energy e is integrated to give a single 
value measure. Currently, the integration is over the full 
horizontal domain and up to approximately 125 hPa. 

The moist static energy is obtained by choosing e 
0.0. However, choosing e = 1.0 gives a relatively large 
weighting to the q component of the norm. In fact with 
this choice the q component of the total energy norm is 
considerably larger than other components. A number 
of challenges remain in terms of properly utilizing moist 
observations and representing moisture in the linear 
model. It may not be sensible to give moisture its full 
weight in the metric before properly investigating the 
behavior. Doing so may lead to misinterpretation of the 
observation impacts, for example by heavily skewing 
impacts toward moist-sensitive instruments. However, 
since the choice of metric is largely arbitrary, it is pos- 
sible to adjust e and tune the relative weighting of the q 
component. Experiments presented here that employ 
the moist norm use a value of e = 0.3; this produces 
approximately equal weighting between the tempera- 
ture and specific humidity components of the norm. 
Experiments using the dry total energy use e = 0. 

f. Filtering of problematic profiles 

Figure 1 shows an adjoint integration initialized with 
the dry total energy at 0000 UTC 18 March 2012. Lin- 
earized moist physics is switched on in the model. The 
adjoint is propagated 6h to 1800 UTC 17 March 2012. 
Figure 1 shows the sensitivity with respect to specific 
humidity at this time (i.e., the end of the adjoint run). 
For the most part the model is well behaved and sensi- 
tivities to moisture are seen in reasonable places, espe- 
cially where storms or fronts are occurring. However, 
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a large localized sensitivity is seen over northwestern 
North America. Sensitivity in the tropics is present, al- 
though it is masked by the dominant behavior in the 
plot. Further integration of the adjoint causes this domi- 
nant sensitivity to grow to a very large magnitude. 

When developing a linear scheme, it is important to 
consider the stability. There are two issues regarding 
stability: first, as to whether the scheme is paired with 
a proper time-stepping scheme with a small enough time 
step and, second, as to whether the scheme contains 
growing modes. Generally, if the linear scheme is not 
paired with a suitable time step and temporal discretiza- 
tion, it will be apparent very quickly in a diverging solu- 
tion. Further, it is generally quite rare for the nonlinear 
scheme to be temporally stable and the linear version of 
the scheme to not be. On the other hand, the presence of 
growing modes can often be an issue in the linearized 
version of the model. For example, consider baroclinic 
instability in the atmosphere, represented by growing 
modes in the nonlinear system. Later on, equilibrium is 
returned in the atmosphere by other processes, perhaps 
represented in a separate part of the nonlinear model. 
A simplified linear model that can represent the per- 
turbations resulting from the baroclinic instability but 
that does not represent the processes responsible for 
the restoration to equilibrium would quickly run into 
trouble. In addition to this the linearization can result 
in large growth rates, for example when a division by 
something small is introduced through the differentiation. 

The large gradients that are observed here appear 
suddenly during just one time step, can first appear at 
various times, and do not seem to be affected by the 
choice of time step; both 15 and 20 min have been tested. 
It seems likely that they are therefore due to the mod- 
eling of instability or large growth introduced in the 
linearization, rather than an issue with the choice of time 
step or time-stepping scheme. 

A way to check for large growth is to examine the 
operator of the tangent linear model, as performed by 
Holdaway and Errico (2013). That study revealed the 
linearized RAS scheme to be stable or close to stable 
everywhere; however, the study only examined two spe- 
cific times. If a large growth is encountered at just one 
particular time and location, it can remain part of the 
solution and lead to further growth, likely seen here. 
Further, even the rather small growth rates seen by 
Holdaway and Errico (2013) could amalgamate over 
time and lead to an issue. 

It is not uncommon to encounter problematic profiles 
such as these when linearizing convection (Errico and 
Raeder 1999; Lopez and Moreau 2005). Generally, some 
kind of filtering must be developed, for which there are 
two general approaches: 


• identify the cause and adapt the linearized (or, better, 
nonlinear) model to prevent problems occurring and 

• use the trajectory to diagnose when the problem 
occurs and ignore or reduce the perturbation there. 

Each of these approaches comes with advantages and 
disadvantages. Making a general change to the linear 
model may be more numerically efficient but can be 
difficult to implement. If the onset of problems is not 
particularly sudden, it can be difficult to identify why 
they occur and therefore produce a targeted enough 
filtering. This could result in too widespread of a cor- 
rection, reducing the closeness to the actual perturba- 
tion trajectory. Identifying problem profiles as they 
occur can more easily produce a focused filtering but 
will be less numerically efficient since it will increase 
the number of calculations. 

Here, the approach is to diagnose problematic profiles 
using the trajectory. Diagnosis must be done using the 
trajectory, rather than the perturbations themselves, so 
as not to cause a discrepancy between the tangent linear 
and adjoint models. 

The tangent linear model is written as 

y' = Mx\ (ll) 

where the vectors y' and x' represent the perturbation 
variables at the end and start times, respectively. The 
matrix M is the tangent linear model operator matrix 
and effectively gives the sensitivity of the scheme with 
respect to the input variables; M depends only on the 
trajectory, or reference, variables. The operator of the 
adjoint model is the transpose, M T . 

If large values suddenly appear in the perturbation 
quantities, the implication would be that the linear op- 
erator matrix M contains a large element for that profile. 
Therefore, a problematic location could be filtered by 
computing the linear operator matrix from the tangent 
linear or nonlinear model and identifying an unusually 
large element, or eigenvalue, as in Errico and Raeder 
(1999). However, computing M for every convective pro- 
file would be prohibitively expensive to do. Computing 
a column of M would generally require either an inte- 
gration of the nonlinear model or an integration of the 
tangent linear model. Computing eigenvalues would be 
very demanding. 

Holdaway and Errico (2013) showed that the struc- 
ture of M is relatively simple and that the locations of 
its dominant features can be understood in terms of 
the profiles of temperature and moisture. Based on those 
findings, it should be possible to target specific columns 
of M and filter based on the values of the gradient in just 
those columns. 
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Fig. 2. The typical structure of the forward operator of the linear RAS scheme: (a) SHIdd, (b) SHIdq, (c) dQldt), and (d) SQIdq. The black 
lines enclose the column of the operator that would be calculated for this particular profile during the filtering. The profile shown is one of 
deep convection. Comparisons of the maximum value in the linear heating and moistening rates throughout the layer for every profile 
against the maximum value in the columns of the operator that are computed for that profile: (e) H' vs the column of 3 Hldd, (f) H' vs the 
column of SHISq, (g) Q' vs the column of SQ/SS , and (h) Q' vs the column of SQIdq. Red points show the profiles that would be adjusted if 
the top 2% are filtered, orange and red if 5% are adjusted, and green, orange, and red if 10% are adjusted. Blue points are all the remaining 
profiles. 


Consider only the potential temperature and specific 
humidity variables. For that case the matrix M in Eq. 
(11) is 

/dH dH\ (r) 

~36 Jq 

M = , (12) 

dQ dQ 

\ dd dq ) 

where H = dOldt and Q = dq/dt are, respectively, the 
heating and moistening rates produced by the nonlinear 
scheme and t is the model time step. The full M would 
also contain components corresponding to the wind 
speeds and surface pressure. For the discrete model 
each component of M is itself a matrix and has di- 
mension n X n , where n is the number of model levels. 

Figures 2a-d show the four components of M: (i) 
dH/dO , (ii) dHIdq, (iii) dQ/dO , and (iv) dQ/dq for a par- 
ticular atmospheric profile. The plotted operator is 
equivalent to Fig. 2 in Holdaway and Errico (2013), 
expect here it is computed using the exact tangent linear 
model. Each column of M is successively computed by 
initializing the tangent linear model with a vector of 
inputs x that is zero everywhere, except at the level 
corresponding to the column being computed, where it 
is set to one. The operator corresponds to a profile that 
exhibits deep convection. It is less efficient to compute M 


using the tangent linear model, rather than the nonlinear 
model, but it avoids the possibility of a switch in the nu- 
merics masking a problem and is found to produce a more 
reliable filtering. 

Columns of M correspond to the level and perturba- 
tion variable being multiplied by that column. So if M 
has structure in a specific column, it describes a sensi- 
tivity to perturbations in that variable at that level. The 
rows of the column correspond to the levels at which 
the response to the perturbation occurs. 

Figures 2a and 2c give the heating rate and moistening 
rate sensitivity with respect to temperature, respec- 
tively. For this profile the dominant sensitivity for tem- 
perature is to perturbations at levels 70 and 45 in the 
model, noted by the structure in those columns in Figs. 
2a and 2c. Level 70 is the level of the cloud base and level 
45 is where the nonlinear heating rate and upward mass 
flux are at their maximum. For the heating rate the 
response to perturbations at these two levels occurs 
throughout the convective region, seen by the struc- 
ture across rows 45-62 in Fig. 2a. For the moistening 
rate the response is dominant at level 62, seen by the 
structure in this row in Fig. 2c. Level 62 is the location 
where the moistening rate is maximum. The diagonal 
feature in Figs. 2a and 2d represents a sensitivity to the 
calculation of dry and moist static energies (Holdaway 
and Errico 2013). The positive gradients below and 
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negative gradients above show how a perturbation at 
a given level causes an increase in heating and moist- 
ening at the level below and a decrease at the level 
above, changing the upward transport of temperature 
and moisture. Figures 2b and 2d give the sensitivity 
with respect to specific humidity. For moisture the 
dominant sensitivity is to perturbations in the subcloud 
layer, denoted by the structure in the columns to the 
right in these panels, again the response occurs through- 
out the convective region in the heating rate and at level 
62 in the moistening rate. See Holdaway and Errico 
(2013) for a complete physical interpretation of these 
structures. The kind of behavior seen for this profile is 
evident for a wide variety of profiles. 

The black lines in Figs. 2a and 2c show the column of 
M that corresponds to the level where the nonlinear 
heating rate is maximum. The black lines in Figs. 2b 
and 2d show the column of M corresponding to specific 
humidity one level below the cloud-base layer, inside the 
subcloud layer. Obtaining the sensitivity to perturba- 
tions of specific humidity somewhere in the subcloud 
layer and perturbations of temperature where the heat- 
ing rate is maximum (i.e., the structure enclosed by the 
black lines) will give a large portion of the important 
dominant features. Two integrations of the tangent lin- 
ear model will not be too computationally demanding. 
The filtering is performed by considering these columns 
and whether the values within them are of a reasonable 
size. 

To know whether filtering based on just the two col- 
umns of M highlighted by the black lines in Figs. 2a-d 
will be successful, the magnitudes are examined against 
the linear growth rates for each convective profile si- 
multaneously when running the tangent linear model. 
The tangent linear model is initialized using an analysis 
increment, which is indicative of the kinds of pertur- 
bation magnitudes that can be expected in practice. In 
Figs. 2e-h the maximum value in the computed col- 
umns of M are scattered against the maximum value in 
the perturbation heating and moistening rates: H' = 
dO'/dt and Q! = dq’ldt. For one time step Fig. 2e shows 
the maximum absolute H ', irrespective of which level 
the maximum occurs at, versus the maximum in the 
column of BH/Bd corresponding to the level where 
the nonlinear heating rate H is at its maximum (e.g., the 
column highlighted in Fig. 2a). Figure 2f shows H' ver- 
sus BHlBq, Fig. 2g shows Q' versus BQ/B9 , and Fig. 2h 
shows Q' versus BQ/Bq. It is clear from these figures that 
the maximum value in the two columns of the operator 
matrix increases as the maximum value of the pertur- 
bation heating and moistening rates increases. This 
positive correlation means that filtering profiles based 
on the operator matrix, obtained with just these two 


perturbations, should be possible. In the figures the 
different colors show the profiles that are associated 
with the largest 2% (red), 5% (orange), and 10% (green) 
of elements in the reduced operator matrix. 

To perform the filtering, four constants, for each 
quadrant of the operator matrix, are determined. If 
the maximum value of the computed columns is larger 
than any of the corresponding constants, then that profile 
is filtered. The four constants are chosen by examining 
the computed columns for a 24-h period (72 time steps) 
and then remain fixed for subsequent experiments. 
Values are chosen so as to filter on average around 4% 
of profiles per time step. This ensures all of the prob- 
lematic profiles are dealt with while minimizing the num- 
ber of profiles that are altered. In Figs. 2e-h approximately 
all of the red points and some of the orange points are 
filtered. 

The filtering targets profiles for which H' is largest. 
Therefore, some profiles where strong convection is 
occurring will not be included. Since these will likely be 
important locations, rather than just assuming that no 
convection is occurring, the perturbation quantity is 
reduced by a factor of 10. This retains the sign of the 
perturbation while also preserving some sensitivity to 
the convection that is occurring. Fortunately, since only 
around 4% of profiles need to be filtered, much of the 
convective behavior is not affected by filtering. 

The linearized moist physics has been tested with a 1° 
horizontal resolution (*“110 km at the equator) and with 
72 levels in the vertical. Plans are under way to increase 
the horizontal resolution of the linear model to V 2 0 ; it is 
possible that the amount of filtering required will de- 
pend on the resolution and this will be tested. The linear 
model with moist physics has been tested with 20- and 
15-min time steps. In both cases problematic profiles 
were encountered and had to be filtered. The time step 
did not impact the amount of filtering that was required 
in order to obtain a satisfactory solution. Different values 
of e have been tested in the energy norm. Again, different 
choices did not impact the amount of filtering required. 

3. Validation 

The linear model is validated by considering how 
well it captures the nonlinear perturbation trajectory. 
The nonlinear perturbation trajectory is obtained by 
taking the end-time difference between two integrations 
of the nonlinear model. In one case the initial conditions 
are perturbed (by Ax) and in one case they are not. The 
same perturbation is then used as the initial conditions 
of the tangent linear model to obtain the linear pertur- 
bation trajectory. By definition of the tangent linear 
model, as the size of the perturbation is reduced, the 
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Fig. 3. The nonlinear perturbation trajectory for the (a) virtual temperature T v and (b) specific humidity q at 500 hPa (level 50 in the 
model) and after a 6-h integration beginning at 0000 UTC 17 Mar. (c),(d) The difference between the nonlinear perturbation trajectory 
and the tangent linear model perturbation trajectory after 6h with moist physics switched off. (e),(f) The difference when the moist 
physics is switched on; (left) virtual temperature and (right) specific humidity. 


nonlinear and linear perturbation trajectories should 
converge. Mathematically this is expressed as 


lim ■"(» + - mW . 1 

Ax— »o MAx 


(13) 


where m represents the nonlinear model and x is the 
nonlinear model variables. The numerator gives the 
nonlinear perturbation trajectory and the denominator 
the linear perturbation trajectory. 

In practice Eq. (13) will not hold, even for very small 
initial perturbations, due to switches and nonlinearity 
in the model. However, it is the behavior of the linear 
model in the presence of the realistic perturbation Ax 
that is of principle interest for almost all tangent linear 
and adjoint applications. The kind of perturbation that 
will be encountered in realistic applications is obtained 


by using an analysis increment (i.e., the analysis minus 
the background, Ax = x fl — x fo ). 

For a given perturbation, the error in the linear per- 
turbation trajectory is 

TLM g = MAx — [m(x + Ax) — m(x)] . (14) 

If Eq. (14) is zero, then the tangent linear model captures 
all the details of the nonlinear perturbation trajectory. 

Figures 3a and 3b show the nonlinear perturbation 
trajectory for virtual temperature and specific humid- 
ity at 500 hPa (model level 50) after a 6-h integration. 
The integration is initialized at 0000 UTC 17 March 
2012. Figures 3c and 3d show TLM e , the tangent linear 
model compared to the nonlinear perturbation differ- 
ence, for virtual temperature and specific humidity 
with the moist physics switched off in the linear model. 
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Fig. 4. As in Fig. 3, but for a 24-h integration and only showing specific humidity differences. The difference when the moist physics is 

switched (left) off and (right) on. 


Figures 3e and 3f show TLM e with the moist physics 
switched on. Contours occurring in the difference plots 
correspond to where the structure is not captured by 
the tangent linear model. 

Comparing Figs. 3d and 3f, it is clear that the inclusion 
of moist physics in the linear model significantly reduces 
the difference between the nonlinear and linear specific 
humidity perturbation trajectories. There are a number 
of regions, especially around the tropics, where the mag- 
nitude of the error has been reduced. The largest dif- 
ferences are seen over the Pacific, Indian, and Atlantic 
Oceans and over Southeast Asia. Regions where the 
dry model fails to capture aspects of the perturbation 
trajectory are improved once moist physics are included. 
Other levels are compared systematically by computing 
correlation coefficients between the nonlinear and linear 
perturbation trajectories. At every level below 100 hPa 
the correlation in the specific humidity field was im- 
proved; at levels above, the change is negligible. For 
the levels from 100 hPa to the surface the correlation 
between the nonlinear and linear perturbation trajec- 
tories increases from an average of 0.66 to 0.73 when 
going from dry to moist. At levels at around 100 hPa, 
where deep convection is dominant, the correlation in- 
creases by as much 35%. 

Although there is significant improvement in the cor- 
relation of the specific humidity perturbation trajectory 
at 6 h, there is less improvement or change to the vir- 
tual temperature perturbation trajectory (Fig. 3c ver- 
sus Fig. 3e). The dry-physics configuration produces a 
better representation of the temperature field pertur- 
bation than it does for moisture; correlations for the 
dry tangent linear model average around 0.73 for below 
100 hPa. There are some improvements when switching 
on moist physics, notably off the coast of South Africa, 
over the western Pacific, and over the Southern Ocean. 
Despite seeing little change in the features, the correla- 
tion coefficient is improved at almost all model levels 
in the temperature field. For the moist linear model 
the average correlation is around 0.76 below 100 hPa. 


There is also a minor improvement in correlation for 
the wind fields with moist physics included. Positive 
impact on the wind fields results from the linear mod- 
eling of cumulus friction. 

That the linearized moist physics model performs well 
for the 6-h window suggests it will prove useful in vari- 
ational data assimilation applications. In the current 
configuration used at GMAO a 6-h window is used. In 
the proposed 4DVAR system the linear perturbation 
trajectory will be used across this window. Further, the 
observation operator employs the linear model to pro- 
duce model space equivalents to the observations; ac- 
curate representation of moisture in the linear model 
over the analysis window is essential for assimilating 
moisture-affected observations. 

Figure 4 shows TLiVl,, at 500 hPa for a 24-h integration. 
As for the 6-h integration the inclusion of linearized 
moist physics results in the largest difference in the mois- 
ture field (temperature not shown). The most significant 
improvement is seen over the Pacific Ocean, over South- 
east Asia, and off the southeast coast of Africa. 

For the 24-h integration the correlations are much 
lower than they are for 6 h. For the dry model the aver- 
age temperature correlation below 100 hPa is around 
0.32 and the average moisture correlation is around 
0.17. When including the linearized moist physics, the 
temperature correlation improvement is neutral and 
the moisture correlation increases to 0.2. These lower 
correlations are expected since the nonlinearities will 
cause a drift from the nonlinear perturbation trajectory. 
However, there is some improvement at almost every 
model level, either for the temperature or the moisture. 
For the wind fields the correlations are largely un- 
changed over 24 h. 

The maximum perturbation over 24 h occurs for the 
moist model over Western Australia, seen in Fig. 4b. 
The magnitude is a little larger than expected and likely 
results from instability that is not captured by the filter- 
ing. However, it is clear that the overall difference be- 
tween the linear and nonlinear perturbation trajectories 
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is decreased when including moist physics and that 
problematic points are being successfully filtered without 
the overall solution being disrupted. If filtering is switched 
off (not shown), a number of locations have very un- 
realistically steep gradients. 

Having an occasional, isolated, slightly too large value 
should not have a negative impact on the applications 
of the tangent linear model. For example it would not 
modify the observation impacts. For the 6-h integration 
the largest values produced by the moist model are sim- 
ilar to the largest values produced by dry model, sug- 
gesting the filtering is sufficient. 

The dot product test is used to ensure that the adjoint 
and tangent linear models have been coded in an equiva- 
lent way. This is done by checking that y T (Mx) = (M T y)x. 
The same level of similarity that is encountered for the 
dry model is found with the moist physics switched on. 

4. Observation impacts 

The GMAO routinely computes adjoint-based ob- 
servation impacts to monitor the large network of ob- 
servations and instruments (Gelaro et al. 2010). The 
adjoint-based observation impact tool is based on the 
work of Langland and Baker (2004). Employing adjoint- 
based impacts is a very powerful and useful tool. Im- 
pacts can be examined per instrument, per channel, for 
different regions of the globe, in a time series, and in 
averages. Metrics that are available include impact per 
analysis, impact per observation, fraction of beneficial 
observations, and observation count per analysis. 

Impacts are computed by integrating two free-running 
forecasts over a 24-h window. One forecast is initialized 
using the analysis x fl and one is initialized using the 
background x b . The forecast initialized with the anal- 
ysis benefits from an extra set of observations and so 
will have a smaller error at the end time. The error at 
the end time is given by 

e f = (/ - x') T P T EP(xf - x r ) , (15) 

where E is a matrix that defines the energy norm [using 
Eq. (10)] and P is used to select the domain over which 
the error is calculated. Superscript / denotes the fore- 
cast and superscript t denotes the truth, or verification. 
The truth is approximated from the model analysis at 
the verification time. The nonlinear observation impact is 
then given by the difference between the errors for the 
two forecasts: e/(x^) — e/(x^). It is the reduction in error 
due to the extra observations and analysis. 

The adjoint is used to propagate the energy norm 
gradient backward 24 h and obtain sensitivities at the 
beginning of the window. The sensitivities are passed 


through the adjoint of the data assimilation system to 
convert them to observation space and the impacts. 
The algorithm for estimating the observation impact, 
e /( x a) ~ e f( x i , )- us ' n 8 the linear model is described in the 
appendix in Langland and Baker (2004). In observation 
space the impact estimate is given by the vector product: 


^(y ” Hx fe ), k t 


( d A+ d A) 

\ dX a dX b) 


(16) 


where y are the observations, H is the linearized obser- 
vation operator, and K T is the adjoint of the data as- 
similation procedure. Sensitivities djf/dx a and dJ b ldx b 
are the gradients of cost functions describing the error 
in the two forecasts, then mapped to observation time 
using the adjoint. 

How good of an approximation the linear observation 
impacts give of the full nonlinear observation impact 
depends on how good of an approximation is obtained 
of djf/dx a and dJ b /dx b and on any approximation made 
in the data assimilation system. There are a number of 
key factors involved. First, how much of the model is 
represented in the linear model: missing physics will 
diminish the approximation? Second, how linear is the 
nonlinear model? (Nonlinearity will cause the linear 
model to drift from the perturbation trajectory.) Third, 
how many approximations are made in the methodology 
and are aspects of the methodology accounted for in 
the adjoint model? For example, the GEOS-5 data as- 
similation uses incremental analysis updating (Bloom 
et al. 1996), meaning observations are applied to the state 
gradually over the assimilation window. The adjoint does 
not take this into account, reducing its accuracy. 

As described by Gelaro et al. (2010), the NASA 
GEOS-5 model uses a double outer loop in the atmo- 
spheric data assimilation. Tremolet (2008) showed that 
systems utilizing multiple outer loops require the use of 
a second-order adjoint model in order to properly cap- 
ture the observation impact. A second-order adjoint is 
not used in the GEOS-5 data assimilation system, so this 
would present a potentially large approximation in the 
methodology. To circumvent this, the assimilation is 
performed using a single outer loop and the same min- 
imization algorithms are used for both the analysis 
(forward) and the sensitivity (backward) parts of the 
integration. The observation impacts computed using 
this single outer loop mode are compared with the op- 
erational observation impacts and are found to be in 
good agreement. 

Figure 5 shows the global nonlinear and linear ob- 
servation impacts for a month, from 17 March until 
17 April 2012. In Fig. 5 the two positive curves show the 
forecast errors ef for the forecasts initialized from the 
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Fig. 5. The forecast error measured in the dry energy norm. The 
positive solid curve shows the error for forecasts initialized using 
the analysis, and the positive dashed curve shows error for forecasts 
initialized using the background. The negative dashed curve shows 
the difference between the forecast errors: the nonlinear obser- 
vation impact. The negative gray curve shows the total adjoint 
impact when using the dry model; the negative black curve shows 
the adjoint impact when moist physics are included. 

analysis and background states. The forecast error is 
computed once daily at the 0000 UTC time. Errors are 
computed using the standard dry energy norm. The 
dashed negative curve shows the nonlinear observation 
impact and is the difference between the two positive 
curves. The two solid negative curves show the linear 
observation impacts. The gray curve shows the linear 
observation impact when using only the dry physics in 
the adjoint, and the black curve shows the linear ob- 
servation impact when including moist physics in the 
adjoint. 

With moist physics included in the adjoint the amount 
of the nonlinear observation impact captured by the 
linear model increases from approximately 77.54% on 
average to 82.59%. This difference between the dry 
model and the moist model is statistically significant 
at the 95% confidence interval. 

Figure 6 shows the average total impact per analysis, 
for each instrument currently assimilated at the GMAO. 
The average total impact per analysis shows how much, 
on average across the month-long period, each instru- 
ment contributed to the total observation impact, shown 
in Fig. 5. The sum of all the blue bars gives the average 
value of the negative black curve in Fig. 5; the sum of 
the red bars is the average of the negative gray curve. 
The observation impacts are computed for the same 
time period (17 March-17 April 2012). The figure com- 
pares the impact for each instrument when using the dry 
adjoint model (red bars) and when moist physics are 
included in the adjoint (blue bars). In both cases the 
error is measured using the dry energy norm so as to 
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Fig. 6. The 24-h forecast observation impacts per analysis for 
each instrument for the period 17 Mar-17 Apr 2012. Red bars show 
the dry-physics configuration of the adjoint and the blue bars show 
the moist-physics configuration. A dry norm is used. 

estimate only the effect of adding the linearized moist 
physics. 

Table 1 numerically compares the impacts found using 
the dry physics and those found with the moist physics. 
It shows the percentage of the total impacts using moist 
physics that is found when using dry physics (i.e., 100 X 
red bars/blue bars in Fig. 6). Values larger than 100% in 
Table 1 represent a reduced reported impact once moist 
physics is used in the adjoint. 

When linearized moist physics are included, all but 
three instruments are reported as having an increased 
positive impact on the forecast. The three instruments 
for which the reported impact is reduced are dropsonde, 
Next Generation Doppler Radar (NEXRAD) winds, and 
profiler winds. These are all instruments that have 
a very small overall impact on the forecast error. The 

Table 1. The percentage of the total observation impact when 
the moist vs dry model is used. The dry norm is used in both cases. 
Percentages in the table effectively show the magnitude of the red 
bars relative to the magnitude of the blue bars in Fig. 6. 


Instrument 

By dry (%) 



Aircraft 

97.43 

MHS 

49.14 

AIRS 

94.70 

MODIS 

37.35 

AMSU-A 

94.03 

NEXRAD 

111.56 

ASCAT 

74.02 

Pilot balloon (pibal) 

87.90 

Dropsonde 

136.42 

Profiler wind 

111.14 

GPS radio occultation 

86.23 

Radiosonde 

95.79 

(GPSRO) 

HIRS 

90.53 

Satellite wind 

93.15 

IASI 

97.59 

TMI rain 

42.80 

Land surface 

75.50 

WindSat 

78.78 

Marine surface 

86.58 

Total 

94.09 
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Fig. 7. The total impact per channel for the AIRS instrument. 
The time period is the same as in Fig. 6. The impacts are generated 
using dry physics and a dry norm. 


Fig. 8. As in Fig. 7, but computed using moist physics in 
the adjoint. 


six largest overall impacts come from aircraft, the Ad- 
vanced Microwave Sounding Unit-A (AMSU-A), the 
Atmospheric Infrared Sounder (AIRS), the Infrared 
Atmospheric Sounding Interferometer (IASI), radio- 
sondes, and satellite winds. For all of these instruments 
the reported impact is increased when moist physics is 
included; the dry configuration captures between 91.15% 
and 97.59% of the impact. The instruments that have 
the largest differences between the dry and moist configu- 
rations are Advanced Scatterometer (ASCAT) winds, land 
surface (land stations from the Meteorological Assimila- 
tion Data Ingest System), Microwave Humidity Sounder 
(MHS), Moderate Resolution Imaging Spectroradiometer 
(MODIS) winds, and Tropical Rainfall Measuring Mission 
(TRMM) Microwave Imager (TMI) rain rate. For those 
instruments the configuration with dry physics in the ad- 
joint captures between 37% and 75% of the impact. These 
are instruments that directly measure moist processes. 
Instruments such as AIRS and HIRS have moist-sensitive 
channels. However, the channels sensitive to temperature 
give the dominant impact so that the dry model captures 
94.70% and 90.53% of the overall impact. 

Figures 7 and 8 show the total impact per channel for 
the AIRS instrument. Figure 7 shows the impact when 
using the dry-physics model and Fig. 8 shows the impact 
when using the moist-physics model. 

AIRS is an infrared sensor that measures over 2378 
spectral channels, of which approximately 120 are as- 
similated. The channels that are most sensitive to mois- 
ture are between 160 and 190 in Figs. 7 and 8. When the 
moist configuration is implemented, the impact being 
reported from these channels is increased, as would be 
expected. However, the impact is typically 40% as 
large as the impact of the temperature-sensitive chan- 
nels. The temperature-sensitive channels also have an 


increased impact when the moist physics are included 
in the model. 

The HIRS and MHS instruments were also examined 
by channel (not shown). Again, the channels associated 
with measuring humidity in the atmosphere were found 
to have a larger impact when the moist configuration is 
implemented. 

Moist norm case 

So far, only the dry energy norm has been considered 
when computing the impacts. It is of interest to examine 
how the use of a moist norm would affect the observation 
impacts. 

The above experiment is repeated for four config- 
urations in total: dry adjoint model initialized with 
a dry norm (red bars in Fig. 6), dry model moist norm, 
moist model dry norm (blue bars in Fig. 6), and moist 
model moist norm. Table 2 shows the percentage of 
the nonlinear error captured by the four linear model 
configurations. 

It is evident from Table 2 that including moist physics 
in the adjoint model increases the fraction of the ob- 
servation impact captured in the linear approximation. 
This suggests that a better approximation of d/{/dx fl 
and djfydxt, in Eq. (16) is obtained. For the dry norm 
(e = 0.0) case the percentage captured by the linear 
model increases by around 5% when including moist 


Table 2. Comparison of the percentage of the nonlinear impact 
captured by the linear configurations. 



Dry norm (%) 

Moist norm (%) 

Dry-physics model 

77.54 

77.33 

Moist-physics model 

82.59 

80.01 


428 


MONTHLY WEATHER REVIEW 


Volume 142 


WINDSAT Wind 
TMI Rain Rate 
Satellite Wind 
Radiosonde 
Profiler Wind 
PIBAL 
NEXRAD Wind 
MODIS Wind 
MHS 

Marine-Surface 
Land-Surface 
IASI 
HIRS 
GPSRO 
Dropsonde 
ASCAT Wind 
Aqua AIRS 
AMSUA 
Aircraft 

- 0.5 - 0.4 - 0.3 - 0.2 - 0.1 0.0 

Observation Impact (J/kg) 

Fig. 9. As in Fig. 6, but comparing the dry-physics configuration 
with a configuration that has moist physics in the adjoint and uses 
the moist norm. 


physics. However, this is a somewhat unfair comparison. 
When including moist physics, a sensitivity will develop 
throughout the adjoint integration in regions where moist 
physics is occurring, even when the initial condition for 
specific humidity is zero. This means moisture is taken 
into account when computing the impact. 

Consider an instrument that measures moisture over 
a location where moist physics is occurring. It is likely 
that a significant portion of its impact will be evident in 
the moisture field and the moist adjoint model can di- 
rectly estimate that impact. Since the dry-physics linear 
model and the dry norm do not directly consider mois- 
ture in their calculations, the impact of that same in- 
strument can only be measured indirectly through the 
other fields. It is therefore not surprising that a large 
increase is observed when comparing the dry and moist 
linear estimates against the nonlinear impact measured 
with a dry norm. 

When the moist norm (e = 0.3) is used, a more modest 
increase of around 3% is seen when going from a dry 
adjoint model to the moist adjoint model. This is a fairer 
comparison since in both cases moisture is taken into 
account in the nonlinear impact calculation and is in- 
cluded in the initial conditions of the adjoint. Even the 
dry linear model has some chance of modeling the im- 
pact on the moisture, through advection of the initial 
specific humidity field. That an increase of around 3% is 
observed here fairly demonstrates that the addition of 
linearized moist physics has increased the accuracy of 
the linear model. 

Figure 9 compares the total impact of the various 
observation systems when using dry physics with the dry 
norm and moist physics with the moist norm for the total 



Total Impact (J/kg) xlCT 2 

Fig. 10. As in Fig. 7, but computed using moist physics in the 
adjoint and using a moist norm. 


impact. Figure 9 shows the percentage captured by the 
dry-dry configuration compared with the moist-moist 
configuration. 

When using the moist physics with the moist norm, the 
reported impacts are much larger than they are when 
using the dry-dry configuration. Of the six instruments 
that produce the largest impact, the dry-dry configura- 
tion only captures between 75% and 90% of the impact. 
The largest change for these instruments is for AIRS, 
for which the dry configuration captures around 75% of 
the moist impact. For the other moist-sensitive instru- 
ments there are also bigger differences. For HIRS around 
66% is captured and for MHS only 17% is captured by 
the dry configuration. 

Figure 10 shows the per channel total impact for AIRS 
when using the moist-physics moist-norm configuration. 
The impact from the humidity measuring channels is 
much larger when moisture is included in the norm. 
Impacts for these channels are of a similar order to 
those of the temperature-sensitive channels. The im- 
pact from the temperature channels in the 100-130 
range is also increased. 

The fourth configuration is the dry model with a moist 
norm. The impacts from that model are similar to those 
found with the moist model and the moist norm. The dry 
model captures 96.78% of the moist-model impact when 
the moist norm is used, compared to 94.09% when the dry 
norm was used. This smaller difference between the dry 
and moist models with the moist norm in terms of im- 
pacts is in agreement with the smaller difference be- 
tween them in terms of the percentage of the nonlinear 
impact captured (Table 2). 

A summary of the percentage of the impact captured 
by the dry model with the dry-norm configuration rela- 
tive to the other configurations is shown in Table 3. 
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Table 3. The total percentage of the impact captured by the dry- 
model dry-norm configuration relative to the other configurations. 


Configuration % 

Moist model dry norm 

94.09 

Dry model moist norm 

85.49 

Moist model moist norm 

82.74 


The two moist-norm configurations report the largest 
overall impacts. By taking moisture into account in 
the norm, the amount of energy that is measured is in- 
creased and this will influence the impacts. An addi- 
tional mechanism by which observations can impact 
the forecast is being measured. For example, an in- 
strument that measures moisture below an intensifying 
storm can have a much larger impact if the moisture is 
represented in the sensitivity, whether it is moderated by 
moist physical processes or not. The inclusion of mois- 
ture in the norm results in many areas where the impact 
due to moisture is being measured. This results in a 
larger relative change to the impacts. 

For all of these experiments a factor of e = 0.3 is in- 
cluded in the moisture term in the norm. If e = 1.0 is 
used instead, then the change in impacts when alter- 
nating between the dry norm and the moist norm is very 
large, especially for moist-sensitive instruments. In the 
current operational model AIRS ranks at about fifth in 
terms of importance for the tropics, after AMSU-A, 
radiodsondes, aircraft reports, and IASI. Even with 
a dry model, using e = 1.0 and the moist norm, AIRS is 
elevated to the most important set of observations here, 
becoming around 25% larger than even AMSU-A. This 
further reinforces the choice to limit the magnitude of 
the moisture term in the energy norm. 

5. Adjoint sensitivity: Storm case study 

An important research-based application of the line- 
arized model is adjoint-based sensitivity analysis. While 
being only a by-product of the observation impact cal- 
culation, the sensitivity fields themselves are of interest 
since they contain information about sensitive regions 
in the initial state. The adjoint initial condition can be 
chosen around an area of interest so that sensitivities 
to just that area are obtained. 

Of particular interest in numerical weather prediction 
is the development of cyclonic systems and the sensi- 
tivity of that development to the initial state. Under- 
standing where the strongest sensitivities lie can be 
helpful in the design of future instruments and in tailor- 
ing observing networks. Most meteorological phenom- 
ena of interest involve both dry and moist dynamics, 
such as hurricanes or midlatitude depressions. Using 


the moist-adjoint model will improve the understanding 
gained from sensitivity studies of these systems. 

As an example of such sensitivities, a strengthening 
depression that is tracking easterly over the North At- 
lantic Ocean toward Europe is identified. The adjoint 
is initialized as the storm reaches its maximum intensity, 
at 0000 UTC 25 March 2012. The energy norm is com- 
puted for a 4° X 4° box around the storm, centered at 
56°N and 40° W. The adjoint is integrated backward 24 h 
to 0000 UTC 24 March 2012. 

Figure 11 shows the adjoint sensitivity at 0000 UTC 
24 March 2012 when using the dry adjoint model and dry 
norm; Fig. 12 shows the adjoint sensitivity when includ- 
ing moist physics in the adjoint and retaining the dry 
norm. Note that there is a change of scale between 
Figs, lid and 12d. The sensitivities are shown for the 
500-hPa level and are shaded red to yellow. Also shown 
in Figs. 11 and 12 is the mean sea level pressure and the 
convective precipitation rate at 0000 UTC 24 March 2012, 
contoured in black and blues, respectively. From 0000 UTC 
24 March to 0000 UTC 25 March 2012 the center of the 
storm moves into the box, also shown in Figs. 11 and 12. 

Comparing Figs. 11 and 12, it is clear that the sensi- 
tivity with respect to specific humidity changes signifi- 
cantly once moist physics are included in the adjoint. 
For the dry model the sensitivity with respect to q is on 
the order of 10 -:> kg kg~ 1 . The specific humidity variable 
itself is generally three orders of magnitude smaller than 
the temperature and wind variables. That the sensitivity 
field with respect to q is of a similar order to the sensi- 
tivity with respect to other fields would mean the over- 
all contribution to the energy would be three orders of 
magnitude smaller than for the other fields, giving the 
impression that total energy is not sensitive to specific 
humidity. When moist physics is included, the sensitivity 
with respect to specific humidity increases by three or- 
ders of magnitude. Given the relative magnitudes of the 
other variables, this would imply that the total energy 
in the box is as sensitive to specific humidity as it is to 
temperature and horizontal wind. 

As well as the relative magnitude, the structure of 
the sensitivity with respect to the specific humidity also 
changes once moist physics is included. Overall, the 
sensitivity is concentrated nearer the low center than it 
is for other fields. There is also less variation with height 
(not shown). At most model levels the sensitivity re- 
mains close to the center of the storm, generally on the 
right-hand side of the storm. The sensitivity also tends to 
lie close to the region where there is convection occur- 
ring, shown by the presence of convective precipitation. 
Trailing behind the storm is a front where strong con- 
vection is evident. The sensitivity is not in the same lo- 
cation as this front for a 24-h adjoint integration. 
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Fig. 11. Adjoint sensitivities to (a) u ' , (b) v' , (c) T' v , and (d) q' at 24 h, shown in red-yellow shading. Sensitivities are shown for ap- 
proximately 500 hPa. Adjoint initialized with energy norm forecast error in a box 54°-58°N, 38.5°-42.5°W, at 0000 UTC 24 Mar 2012. The 
convective precipitation at — 24 h is shown in blue-green contours. In black contours is the sea level pressure at — 24 h. Sensitivities are 
generated using the dry adjoint, i.e., without moist physics turned on. The box shows the region where the forecast error is measured. The 
center of the low pressure system moves into the box over the 24-h window. 


For the sensitivity with respect to it, v, and T m the 
structure of the field is similar with or without moist 
physics included in the adjoint model. The only change is 
a slight increase in the magnitudes when moist physics 
is included. As expected, when higher levels were ex- 
amined, the sensitivity was found to rotate clockwise 
around the storm with increasing height. For high 
enough levels the sensitivity with respect to the wind 
fields stretches out over northern Canada into the large- 
scale synoptic flow features (not shown). 

The sensitivity experiment was repeated for both 
moist and dry norms with moist and dry physics in the 
adjoint. The sensitivities that are obtained when the 
moist-physics model is used are not found to be de- 
pendent on the choice of e. Use of both dry and moist 
norms, with various choices for s, produces approxi- 
mately the same overall sensitivity fields. This is dis- 
cussed further below. 

6. Dependency on moisture in the norm 

When using adjoint techniques, the choice of norm, 
which gives the initial conditions for the adjoint model, 


requires careful consideration. In choosing a specific 
norm a specific question is asked. Choosing the total 
energy norm, for example, poses the question: What is a 
quadratic measure of the perturbation field, that looks 
like energy, sensitive to? Or in the observation impact 
experiments: How much of a role does a particular 
observation play in reducing the error in the system, 
measured in this way? The choice of norm used to ini- 
tialize the adjoint must be tailored to the problem being 
addressed. The energy norm is often used because it 
offers a convenient way of combining different model 
variables that have different units. Using the energy 
norm in the observation impact experiments allows for 
the quantification of the impact that observations have 
on all the variables simultaneously. 

In both the observation impact experiments and the 
sensitivity study, dry and moist energy norms have 
been tested. Clearly, just the choice of e in the norm 
has a large influence on the observation impacts that 
are obtained. However, it is interesting to note that 
for the sensitivity study the counter result is obtained, 
the same sensitivity is found regardless of the choice 
of e. 
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Fig. 12. As in Fig. 11, but with moist physics turned on in the adjoint model. Note the change in scale in (d) compared with Fig. lid. 


In this study it is found that when the moist norm is 
chosen, the modification of the sensitivity with respect 
to q due to the moist-physics schemes dominates over 
the modification due to advection. So a location that 
has moisture sensitivity only due to the initial condi- 
tions, or due to advection of the initial conditions, will 
have a smaller overall sensitivity than a neighboring lo- 
cation that has highly active moist physics. In addition to 
this, it is found that the initial conditions of moisture can 
be quickly forgotten when the linearized moist-physics 
schemes are called. Consider for example Eq. (1); in the 
formulation of q s there is no dependence on specific 
humidity, only temperature and pressure. Much of the 
changes due to the moist physics are computed without 
considering the moisture in the initial conditions. The 
modification of the sensitivity by the moist physics de- 
pends much more strongly on the temperate in the initial 
conditions than the moisture in the initial conditions. 
After a few time steps at a location with active moist 
physics, the schemes will sufficiently modify the sen- 
sitivity to the point that the moisture initial conditions 
provided by the norm are forgotten. 

The difference between the observation impact ex- 
periment and the limited-area sensitivity study is the 
relative number of locations where the moist-physics 
schemes are active. For the sensitivity study the adjoint 


is initialized with an energy norm computed only for the 
region close to the center of the storm. Over the temporal 
and spatial domain being examined both convection 
and large-scale condensation schemes are strongly mod- 
ifying the sensitivity at every location. In addition to this 
the sensitivity quickly spreads out to areas where the 
norm was not calculated, further diluting the effect of the 
initial conditions. 

For the observation impact study the adjoint model is 
initialized with a norm calculated for the whole globe. 
In this global situation there will be many locations 
where the moist-physics schemes are highly active and 
therefore modifying the sensitivity. However, there will 
be far more locations where the norm, and therefore 
the initial conditions, contains moisture sensitivity but 
where moist physics are not active. For these locations, 
where moist physics do not modify the sensitivity, the 
original choice of e will determine the final sensitivity. 
Since there are a large number of locations like this, the 
choice of e has a large overall impact. 

That the limited-area sensitivity study is not sensitive 
to the choice of e is not to say that it is not sensitive to 
the choice of norm entirely. If a different type of norm 
was chosen, for example the root-mean-square error of 
moisture in the box, then a very different sensitivity 
field would likely result. When performing this kind of 
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study, one has to ask what is of interest. It may be that 
a measure like total energy is of most interest or it may 
be some other metric. 

7. Summary 

A linearization of moist physics has been developed 
for use in NASA’s GEOS-5 atmospheric data assimila- 
tion model. Convection is modeled using a linearization 
of the relaxed Arakawa-Schubert scheme. Large-scale 
condensation and precipitation are modeled using a lin- 
earization of a highly simplified scheme that removes 
supersaturation. 

The new linearizations have been validated by com- 
paring the tangent linear perturbation trajectory with 
the nonlinear perturbation trajectory that includes the 
full moist physics. For both 6- and 24-h tangent linear 
integrations the inclusion of moist physics is found to 
improve the correlation between the nonlinear and lin- 
ear perturbation trajectories. In computing the obser- 
vation impacts it is found that including moist physics 
increases the representation of the nonlinear impact by 
the linear model. 

Since including moist physics increases the likeness 
between the linear and nonlinear observation impact, it 
should provide more realism in the impacts. Almost all 
of the instruments assimilated in GEOS-5 are reported 
as having a larger impact when the moist-adjoint con- 
figuration is considered. With the moist physics switched 
on, the estimate of the impact is more accurate since it is 
able to account for a larger fraction of the overall ob- 
servation impact. Moist-sensitive observations from in- 
struments such as MHS, TMI, infrared sounders, and 
some wind observations see the largest changes in the 
impact per analysis. Some instruments see a smaller 
number of their observations being beneficial but the 
difference is small and the total impact can still increase. 
Impacts from the AIRS, HIRS, and MHS instruments 
were examined by channel. Those channels with a sensi- 
tivity to moisture have the largest positive change. Most 
temperature -sensitive channels also see a small positive 
change. 

A case study of an intensifying storm over the North 
Atlantic is used to examine the effect of including moist 
physics on sensitivity studies. With moist physics in- 
cluded, a large sensitivity with respect to the specific 
humidity is encountered. Whereas for the dry model 
the dominant sensitivity is with respect to horizontal 
wind and temperature, for the moist model the sensi- 
tivity with respect to all variables is equally large. 

An important consideration when using the linearized 
adjoint model is the choice of norm. Currently, a dry 
total energy norm is used as the metric. Including moist 


physics in the model raises the question of whether 
moisture should also be considered in the metric. Four 
configurations were considered for both observation 
impacts and the storm sensitivity study. For the obser- 
vation impacts there is a strong dependency on the 
presence of moist physics in the energy norm. Including 
moisture results in a large change in the observation 
impacts, whether the model physics are dry or moist. 
Choosing the natural form of moist static energy results 
in AIRS being reported as the instrument with the 
largest impact, ahead of AMSU-A. To reduce the em- 
phasis on observations affected by moisture, a factor 
(e = 0.3) is included in the moisture term. For the sen- 
sitivity case study the choice of e is not as important. 
Whether a dry or moist norm is chosen, the same overall 
sensitivity field is obtained. For observation impacts 
there are many locations where no moist physics is oc- 
curring. Including moisture in the metric will allow in- 
struments measuring moisture to have an impact at 
those locations. For the sensitivity study the locations 
where moist physics is active are dominant and the ini- 
tial conditions are dominated by the physical behavior. 

It remains to test this suite of linearized moist physics 
within a 4DVAR framework; it would be interesting 
to examine how well observations of moist processes 
can be assimilated with this particular linearized moist- 
physics package. Currently, the large-scale precipitation 
is represented in a very crude manner. It would likely be 
beneficial to try and include more complex processes 
such as reevaporation and to distinguish between types 
of precipitation. 
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